home *** CD-ROM | disk | FTP | other *** search
/ Packard Bell - Internet on a CD / internet on a cd.cdr / Internet / sites / Clementine_NASA / pds1.c < prev    next >
Encoding:
C/C++ Source or Header  |  1998-07-16  |  12.7 KB  |  529 lines
  1. /*
  2.  *      THIS ROUTINE IS PART OF THE CLEMENTINE PDS FILE READER PROGRAM.
  3.  *      IT WAS WRITTEN BY ACT CORP. IN DIRECT SUPPORT TO THE
  4.  *      CLEMENTINE (DSPSE) PROGRAM.
  5.  *
  6.  *      IF YOU FIND A PROBLEM OR MAKE ANY CHANGES TO THIS CODE PLEASE CONTACT
  7.  *      Dr. Erick Malaret at ACT Corp.
  8.  *                      tel: (703) 742-0294
  9.  *                           (703) 683-7431
  10.  *                       email: nrlvax.nrl.navy.mil
  11.  *
  12.  *
  13.  */
  14. #include <stdio.h>
  15. #include <math.h>
  16. /*#include <alloc.h>*/
  17. #include <malloc.h>
  18. /* #include <mem.h> */
  19. #include <string.h>
  20. #include "jpeg_c.h"
  21. #include "pds.h"
  22.  
  23. #define READ_PARAM "rb"
  24. #define WRITE_PARAM "wb"
  25.  
  26. PDSINFO pds;
  27. FILE    *qparm;
  28. char    histfname[80];
  29. int     USEHIST    = 0;
  30. int     OPTIMIZE   = 1;
  31. int     SLOWDECOMP = 0;
  32. extern long     *DCTHist[64];
  33. extern float    *Rn[64];
  34. extern float    Q[64];
  35. extern float    C[64], Ct[64];
  36.  
  37. void init_qt(FILE *fptr);
  38. void readhufftbls(FILE *fptr);
  39. void pds_decomp(FILE *fptr,CHARH *p,long sizej,long sizei);
  40. #ifdef  sun
  41. void PClong2SUNlongVector(unsigned long invec[],int npts) ;
  42. void PCshort2SUNshortVector(unsigned short invec[],int npts) ;
  43. #endif
  44. extern void getDCTHist(BitStream *bs, long rows, long cols);
  45. extern void getRn(void);
  46. extern void initC(float *m, short N);
  47. extern void transp(float *out, float *in, short N);
  48.  
  49.  
  50. PDSINFO *PDSR(char *fname, long *rows, long *cols)
  51. {
  52.     FILE    *fptr;
  53.     int     n;
  54.     int     j,i;
  55.     CHARH   *c;
  56.     char    nstring[84],sdummy[80],buffer[10],low;
  57.     long    sizej;                  /* number of rows in the image */
  58.     long    sizei;                  /* number of columns in the image */
  59.     int     bitpix;                 /* bits per pixel */
  60.     int     record_bytes;
  61.     int     hist_rec,brw_rec,image_rec,ns,brwsize;
  62.     char    cval, *sptr, *ptr;
  63.     char    record_type[20];
  64.     int     COMPRESSED;
  65.     long    k,hdr_size;
  66.  
  67.     if( (fptr = fopen(fname,READ_PARAM)) == NULL){  /* open disk file */
  68.         printf("Can't open %s.\n",fname);
  69.         return NULL;
  70.     }
  71.  
  72.     /* initialize some basic variables */
  73.     bitpix                  = 0;
  74.     sizej=sizei             = 0;
  75.     pds.browse_nrows        = 0;
  76.     pds.browse_ncols        = 0;
  77.     pds.image_nrows = 0;
  78.     pds.image_ncols = 0;
  79.     hist_rec = brw_rec = image_rec = -1;
  80.  
  81.     /* read header */
  82.     do{
  83.         /* read next line of text */
  84.         for (n=0; (cval=fgetc(fptr)); n++) {
  85.             nstring[n]=cval;
  86.             if(cval=='\n') {
  87.                 if( (cval=fgetc(fptr)) != '\r') ungetc(cval,fptr);
  88.                 nstring[++n]='\0';
  89.                 break;
  90.             }
  91.         }
  92.  
  93.         /* find line's first non-space character */
  94.         for (ns=0; nstring[ns]==' ';ns++);
  95.         sptr = &nstring[ns];
  96.  
  97.         if (strncmp("^IMAGE_HISTOGRAM ",sptr,17)==0) {
  98.             n=sscanf(nstring,"%s = %d", sdummy, &hist_rec);
  99.         }
  100.  
  101.         if (strncmp("^BROWSE_IMAGE ",sptr,14)==0) {
  102.             n=sscanf(nstring,"%s = %d", sdummy, &brw_rec);
  103.         }
  104.  
  105.         if (strncmp("^IMAGE ",sptr,7)==0) {
  106.             n=sscanf(nstring,"%s = %d", sdummy, &image_rec);
  107.         }
  108.  
  109.         if (strncmp("RECORD_TYPE",sptr,9)==0) {
  110.             n=sscanf(nstring,"%s = %s", sdummy, &record_type);
  111.         }
  112.  
  113.         if (strncmp("RECORD_BYTES",sptr,12)==0) {
  114.             n=sscanf(nstring,"%s = %d", sdummy, &record_bytes);
  115.         }
  116.  
  117.         if (strncmp("ENCODING_TYPE",sptr,13)==0) {
  118.             n=sscanf(nstring,"%s = %s", sdummy, buffer);
  119.             if ( strstr(buffer,"JPEG") )
  120.                 COMPRESSED = 1;
  121.             else
  122.                 COMPRESSED = 0;
  123.         }
  124.  
  125.         if (strncmp("LINES ",sptr,6)==0) {
  126.             n=sscanf(nstring,"%s = %d", sdummy, &sizej);
  127.         }
  128.  
  129.         if (strncmp("LINE_SAMPLES",sptr,12)==0) {
  130.             n=sscanf(nstring,"%s = %d", sdummy, &sizei);
  131.         }
  132.  
  133.         if (strncmp("SAMPLE_BITS",sptr,11)==0) {
  134.             n=sscanf(nstring,"%s = %d", sdummy, &bitpix);
  135.         }
  136.  
  137.         if (strncmp("END",sptr,3)==0) {
  138.             if ( *(sptr+3)=='\n' || *(sptr+3)==' ' || *(sptr+3)=='\r' ) {
  139.                 hdr_size = ftell(fptr);
  140.                 break;
  141.             }
  142.         }
  143.     } while(1);
  144.  
  145.     /**************   read histogram  ***************/
  146.     if ( hist_rec != -1 ) {
  147.         fseek(fptr,hist_rec-1,0);
  148.         pds.hist = (long *)malloc(256*sizeof(long));
  149.         if(pds.hist == NULL) {
  150.             printf(" histogram memory not allocated \n");
  151.         }
  152.         if ( pds.hist ) {
  153.             fread(pds.hist,sizeof(long),256,fptr);
  154. #ifdef sun
  155.             PClong2SUNlongVector((unsigned long*)pds.hist,256);
  156. #endif
  157.         }
  158.     }
  159.  
  160.     /**************   read browse image **********/
  161.     if ( brw_rec != -1 ) {
  162.         pds.browse_ncols        = sizei/8;
  163.         pds.browse_nrows        = sizej/8;
  164.         fseek(fptr,brw_rec-1,0);
  165.         brwsize = (sizej/8) * (sizei/8);
  166.         pds.brw_imag = (unsigned char *)malloc(brwsize);
  167.         if ( pds.brw_imag )
  168.             fread(pds.brw_imag,sizeof(char),brwsize,fptr);
  169.     }
  170.  
  171.     /*************   read image data ***************/
  172.     if (strncmp(record_type,"UNDEFINED",9)==0) {
  173.         record_bytes=1;
  174.         fseek(fptr,(image_rec-1),0);
  175.     } else {
  176.         fseek(fptr,(image_rec-1)*record_bytes,0);
  177.     }
  178.  
  179.     switch (bitpix) {
  180.     case 8:
  181.         c = (CHARH *)MALLOC(sizej*sizei);
  182.         if ( c == NULL ) {
  183.             printf("Can't allocate memory for image array.\n");
  184.             fclose(fptr);
  185.             return NULL;
  186.         }
  187.  
  188.         if ( COMPRESSED ) {
  189.             qparm = fopen("paramtrs.dat","w");
  190.             init_qt(fptr);
  191.             readhufftbls(fptr);
  192.             pds_decomp(fptr,c,sizej,sizei);
  193.             fclose(qparm);
  194.         } else {
  195.  
  196.             for (j=0, k=0; j<sizej ;j++) {
  197.                 for (i=0; i<sizei; i++) {
  198.                     if ( 1 != fread(&low,sizeof(char),1,fptr) ) {
  199.                         printf("error: possible EOF found at (%d,%d)\n",j,i);
  200.                         break;
  201.                     } else {
  202.                         c[k++] = (unsigned char)low;
  203.                     }
  204.                 }
  205.             }
  206.         }
  207.  
  208.         pds.image = c;
  209.         break;
  210.     default:
  211.         printf("invalid number of bits per pixel\n");
  212.         pds.image = NULL;
  213.     }
  214.  
  215.     /************    Allocate string buffer    **************/
  216.     rewind(fptr);
  217.  
  218.     pds.text = (char *)malloc(hdr_size);
  219.     if ( pds.text ) {
  220.         for (ptr=pds.text,i=0; i<hdr_size; i++) {
  221.             fread(ptr,sizeof(char),1,fptr);
  222.             if ( *ptr == '\r' ) {
  223.                 /* do nothing */
  224.             } else {
  225.                 ptr++;
  226.             }
  227.         }
  228.     }
  229.     *(ptr)='\0';
  230.  
  231.     /*****/
  232.     fclose(fptr);
  233.  
  234.     *rows = pds.image_nrows=sizej;
  235.     *cols = pds.image_ncols=sizei;
  236.  
  237.     return &pds;
  238. }
  239.  
  240.  
  241. /******** Routines that deal with compressed images *******/
  242. float   dfac[8] = {
  243. 0.35355339,0.35355339,0.653281482,0.27059805,
  244. 0.449988111,0.254897789,0.300672443,1.281457724
  245. };
  246.  
  247. void init_qt(FILE *fptr)
  248. {
  249.     short   i;
  250.     unsigned short   scalef,index;
  251.     unsigned short   table[64];
  252.     float   ftable[64];
  253.  
  254.     fread(&scalef,sizeof(short),1,fptr);
  255. #ifdef sun
  256.     PCshort2SUNshortVector((unsigned short*)&scalef,1) ;
  257. #endif
  258.     fprintf(qparm,"tabf: %d\n",scalef);
  259.     fread(table,sizeof(short),64,fptr);
  260. #ifdef sun
  261.     PCshort2SUNshortVector((unsigned short*)table,64);
  262. #endif
  263.     fprintf(qparm,"tabq:\n");
  264.  
  265.     for (i=0; i<64; i++) {
  266.         fprintf(qparm,"%3d ",table[i]);
  267.         if ( (i+1) % 8 == 0 ) fprintf(qparm,"\n");
  268.  
  269.         ftable[i] = (float)(scalef*table[i])/64.0 + 0.5;
  270.         ftable[i] = 4096.0 / (float)floor(ftable[i]);
  271.     }
  272.  
  273.     for (i=0; i<64; i++) Q[i] = ftable[zzseq[i]];
  274.  
  275.     ftable[0] = dfac[0]*dfac[0]*ftable[0];
  276.     ftable[32] = dfac[0]*dfac[1]*ftable[32];
  277.     ftable[16] = dfac[0]*dfac[2]*ftable[16];
  278.     ftable[48] = dfac[0]*dfac[3]*ftable[48];
  279.     ftable[8] = -dfac[0]*dfac[4]*ftable[8];
  280.     ftable[24] = -dfac[0]*dfac[5]*ftable[24];
  281.     ftable[56] = dfac[0]*dfac[6]*ftable[56];
  282.     ftable[40] = -dfac[0]*dfac[7]*ftable[40];
  283.  
  284.     ftable[4] = dfac[1]*dfac[0]*ftable[4];
  285.     ftable[36] = dfac[1]*dfac[1]*ftable[36];
  286.     ftable[20] = dfac[1]*dfac[2]*ftable[20];
  287.     ftable[52] = dfac[1]*dfac[3]*ftable[52];
  288.     ftable[12] = -dfac[1]*dfac[4]*ftable[12];
  289.     ftable[28] = -dfac[1]*dfac[5]*ftable[28];
  290.     ftable[60] = dfac[1]*dfac[6]*ftable[60];
  291.     ftable[44] = -dfac[1]*dfac[7]*ftable[44];
  292.  
  293.     ftable[2] = dfac[2]*dfac[0]*ftable[2];
  294.     ftable[34] = dfac[2]*dfac[1]*ftable[34];
  295.     ftable[18] = dfac[2]*dfac[2]*ftable[18];
  296.     ftable[50] = dfac[2]*dfac[3]*ftable[50];
  297.     ftable[10] = -dfac[2]*dfac[4]*ftable[10];
  298.     ftable[26] = -dfac[2]*dfac[5]*ftable[26];
  299.     ftable[58] = dfac[2]*dfac[6]*ftable[58];
  300.     ftable[42] = -dfac[2]*dfac[7]*ftable[42];
  301.  
  302.     ftable[6] = dfac[3]*dfac[0]*ftable[6];
  303.     ftable[38] = dfac[3]*dfac[1]*ftable[38];
  304.     ftable[22] = dfac[3]*dfac[2]*ftable[22];
  305.     ftable[54] = dfac[3]*dfac[3]*ftable[54];
  306.     ftable[14] = -dfac[3]*dfac[4]*ftable[14];
  307.     ftable[30] = -dfac[3]*dfac[5]*ftable[30];
  308.     ftable[62] = dfac[3]*dfac[6]*ftable[62];
  309.     ftable[46] = -dfac[3]*dfac[7]*ftable[46];
  310.  
  311.     ftable[1] = -dfac[4]*dfac[0]*ftable[1];
  312.     ftable[33] = -dfac[4]*dfac[1]*ftable[33];
  313.     ftable[17] = -dfac[4]*dfac[2]*ftable[17];
  314.     ftable[49] = -dfac[4]*dfac[3]*ftable[49];
  315.     ftable[9] = dfac[4]*dfac[4]*ftable[9];
  316.     ftable[25] = dfac[4]*dfac[5]*ftable[25];
  317.     ftable[57] = -dfac[4]*dfac[6]*ftable[57];
  318.     ftable[41] = dfac[4]*dfac[7]*ftable[41];
  319.  
  320.     ftable[3] = -dfac[5]*dfac[0]*ftable[3];
  321.     ftable[35] = -dfac[5]*dfac[1]*ftable[35];
  322.     ftable[19] = -dfac[5]*dfac[2]*ftable[19];
  323.     ftable[51] = -dfac[5]*dfac[3]*ftable[51];
  324.     ftable[11] = dfac[5]*dfac[4]*ftable[11];
  325.     ftable[27] = dfac[5]*dfac[5]*ftable[27];
  326.     ftable[59] = -dfac[5]*dfac[6]*ftable[59];
  327.     ftable[43] = dfac[5]*dfac[7]*ftable[43];
  328.  
  329.     ftable[7] = dfac[6]*dfac[0]*ftable[7];
  330.     ftable[39] = dfac[6]*dfac[1]*ftable[39];
  331.     ftable[23] = dfac[6]*dfac[2]*ftable[23];
  332.     ftable[55] = dfac[6]*dfac[3]*ftable[55];
  333.     ftable[15] = -dfac[6]*dfac[4]*ftable[15];
  334.     ftable[31] = -dfac[6]*dfac[5]*ftable[31];
  335.     ftable[63] = dfac[6]*dfac[6]*ftable[63];
  336.     ftable[47] = -dfac[6]*dfac[7]*ftable[47];
  337.  
  338.     ftable[5] = -dfac[7]*dfac[0]*ftable[5];
  339.     ftable[37] = -dfac[7]*dfac[1]*ftable[37];
  340.     ftable[21] = -dfac[7]*dfac[2]*ftable[21];
  341.     ftable[53] = -dfac[7]*dfac[3]*ftable[53];
  342.     ftable[13] = dfac[7]*dfac[4]*ftable[13];
  343.     ftable[29] = dfac[7]*dfac[5]*ftable[29];
  344.     ftable[61] = -dfac[7]*dfac[6]*ftable[61];
  345.     ftable[45] = dfac[7]*dfac[7]*ftable[45];
  346.  
  347.     for (i=0; i<64; i++) qtable[i] = ftable[zzseq[i]];
  348. }
  349.  
  350. void readhufftbls(FILE *fptr)
  351. {
  352.     fread(dcbits,sizeof(short),16,fptr);
  353. #ifdef sun
  354.     PCshort2SUNshortVector((unsigned short*)dcbits,16);
  355. #endif
  356.  
  357.     fread(dchuffval,sizeof(char),12,fptr);
  358.     fread(acbits,sizeof(short),16,fptr);
  359. #ifdef sun
  360.     PCshort2SUNshortVector((unsigned short*)acbits,16);
  361. #endif
  362.  
  363.     fread(achuffval,sizeof(char),162,fptr);
  364.  
  365.     inithuffcode();
  366. }
  367.  
  368.  
  369. void pds_decomp(FILE *fptr,CHARH *p,long sizej,long sizei)
  370. {
  371.     BitStream       ibs;
  372.     unsigned short  nb;
  373.     short   i, j, npanels;
  374.     long    nbytes, bytesperpanel;
  375.     short   blocks, rem;
  376.     long    filepos1, filepos2;
  377.     CHARH   *ptr;
  378.     int     FLAG = 0;
  379.  
  380.     cBitStream(&ibs,NULL,INPUT);
  381.  
  382.     filepos1 = ftell(fptr);
  383.     fseek(fptr,0,2);
  384.     filepos2 = ftell(fptr);
  385.     fseek(fptr,filepos1,0);
  386.  
  387.     nbytes = filepos2 - filepos1;
  388.  
  389.     ibs.outstring = (CHARH *)MALLOC(nbytes);
  390.     if ( ibs.outstring ) {
  391.         blocks = 1;
  392.         rem = 0;
  393.         nb = (unsigned short)nbytes;
  394.         if ( nbytes > 60000L ) {
  395.             blocks = nbytes / 32768;
  396.             rem = nbytes % 32768;
  397.             nb = 32768;
  398.         };
  399.         ptr = ibs.outstring;
  400.         for (i=0; i < blocks; i++,ptr+=nb) {
  401.             if ( fread(ptr,sizeof(char),nb,fptr) != nb ) {
  402.                 printf("Error reading data string.\n");
  403.                 FLAG = 1;
  404.                 break;
  405.             }
  406.         }
  407.         if ( rem ) {
  408.             if ( fread(ptr,sizeof(char),rem,fptr) != rem ) {
  409.                 printf("Error reading data string.\n");
  410.                 FLAG = 1;
  411.             }
  412.         }
  413.         ibs.mode = MEMORY;
  414.     } else {
  415.         ibs.bytestream.file = fptr;
  416.         ibs.mode = DISK;
  417.     }
  418.  
  419.     if ( !FLAG ) {
  420.         npanels = sizej/32;
  421.         bytesperpanel = 32*sizei;
  422.  
  423.         if ( OPTIMIZE ) {
  424.             for (i=0; i<64; i++) {
  425.                 DCTHist[i] = (long FAR *)MALLOC(sizeof(long)*513);
  426.                 memset(DCTHist[i],0,sizeof(long)*513);
  427.                 DCTHist[i] += 256;
  428.             }
  429.             for (i=0; i<64; i++) {
  430.                 Rn[i] = (float FAR *)MALLOC(sizeof(float)*513);
  431.                 memset(Rn[i],0,sizeof(float)*513);
  432.                 Rn[i] += 256;
  433.             }
  434.  
  435.             if ( USEHIST ) {
  436.                 FILE *  f1;
  437.                 long    sx, sy;
  438.                 char    cdummy;
  439.                 float   histval;
  440.                 f1 = fopen(histfname,"rb");
  441.                 fread(&sx,sizeof(long),1,f1);
  442.                 fread(&sy,sizeof(long),1,f1);
  443.                 fread(&cdummy,sizeof(char),1,f1);
  444.                 sy -= 256;
  445.                 for (i=0; i<sx; i++) {
  446.                     for (j=-256; j<sy; j++) {
  447.                         fread(&histval,sizeof(float),1,f1);
  448.                         DCTHist[i][j] = (long)histval;
  449.                     }
  450.                 }
  451.                 fclose(f1);
  452.             } else {
  453.                 for (i=0; i<npanels; i++) getDCTHist(&ibs,32,sizei);
  454.                 ibs.BitBuffer = 0;
  455.                 ibs.bytesout = 0;
  456.                 ibs.BitBuffMask = 0x00;
  457.             }
  458.  
  459.             getRn();
  460.         }
  461.  
  462.         for (i=0, ptr=p; i<npanels; i++,ptr+=bytesperpanel) {
  463.             decomp(&ibs,ptr,32,sizei);
  464.         }
  465.  
  466.         if ( OPTIMIZE ) {
  467.             FILE *  outf;
  468.             long    ldummy;
  469.             char    cdummy;
  470.             float   histval;
  471.             outf = fopen("hist.flt","wb");
  472.             if ( outf ) {
  473.                 ldummy = 64;
  474.                 fwrite(&ldummy,sizeof(long),1,outf);
  475.                 ldummy = 513;
  476.                 fwrite(&ldummy,sizeof(long),1,outf);
  477.                 cdummy = 0;
  478.                 fwrite(&cdummy,sizeof(char),1,outf);
  479.             } else {
  480.                 printf("*** Could not open hist.flt ***\n");
  481.             }
  482.             for (i=0; i<64; i++) {
  483.                 DCTHist[i] -= 256;
  484.                 Rn[i] -= 256;
  485.                 if ( outf ) {
  486.                     for (j=0; j<513; j++) {
  487.                         histval = (float)DCTHist[i][j];
  488.                         fwrite(&histval,sizeof(float),1,outf);
  489.                     }
  490.                 }
  491.                 FREE(DCTHist[i]);
  492.                 FREE(Rn[i]);
  493.             }
  494.             if ( outf ) fclose(outf);
  495.         }
  496.     }
  497.  
  498.     if ( ibs.outstring ) FREE( ibs.outstring );
  499. }
  500.  
  501.  
  502. #ifdef sun
  503. void PClong2SUNlongVector(unsigned long invec[],int npts){
  504.     int     i;
  505.     unsigned long   ival,oval;
  506.  
  507.     for(i=0;i<npts;i++){
  508.         ival    = invec[i];
  509.         oval    = ((ival&0x000000ff)<<24) +
  510.               ((ival&0x0000ff00)<<8) +
  511.               ((ival&0x00ff0000)>>8) +
  512.               ((ival&0xff000000)>>24);
  513.         invec[i]= oval;
  514.     }
  515. }
  516.  
  517. void PCshort2SUNshortVector(unsigned short invec[],int npts){
  518.     int     i;
  519.     unsigned short  ival,oval;
  520.  
  521.     for(i=0;i<npts;i++) {
  522.         ival    = invec[i];
  523.         oval    = (ival<<8) + ((ival>>8)&0x00ff);
  524.         invec[i]= oval;
  525.     }
  526. }
  527. #endif
  528.  
  529.